/**
  ETFRoot project, Anyang Normal University && IMP-CAS
  \file yield.cxx
  \brief define cuts to select reaction products
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2022-06-24
  \date 2021-06-24 last modified
  \attention
  changelog
  <table>
  <tr>  <th>Date         <th>Author      <th>Description                   </tr>
  <tr>  <td>2022-06-24   <td>Asia Sun    <td>file created                  </tr>
  </table>

  \copyright Copyright (c) 2021 Anyang Normal U. && IMP-CAS with LGPLv3 LICENSE
*/

const string dir = "/home/asia/Lab/ETFAna/build/";
const string file = "C16TA.root";

const string ff = dir+file;
static const double PI = ETFMath::Pi();

void yield(){
  TFile *f = new TFile(ff.c_str(), "UPDATE");
  TTree *ETF_TREE = (TTree *)f->Get("ETF_TREE");
  if(!ETF_TREE) ETFMsg::Error("main", "cut: nullptr ETF_TREE");
  vector<TCutG *> cutls;

  // CUT A0: choose the incident particle in dE0-tof1 spectrum
	const int n = 100; // number of points for TCutG in the dE1-aoz spectrum
	const char *nucls[] = {"C16", "C15", "C14", "C13", "B15", "B14", "B13", "B12",
    "B11", "Be12", "Be11", "Be10", "Be9", "C16R", "C16U"};
	for(const char *c : nucls) cutls.push_back(new TCutG(c, n));

	for(TCutG *c : cutls){
		c->SetTitle("NUCLIDE CUT");
    c->SetVarY("pid.dE1()");
		c->SetVarX("pid.aoz()");
		c->SetFillStyle(100);
		c->SetLineColor(2);
		c->SetLineWidth(2);
		while(c->GetN()) c->RemovePoint(0);
	} // end for over TCutG objects

	// a: aoz - half length of the ellipse's major axis
	// b: dE1 - half length of the ellipse's minor axis
	vector<double> x0 = {2.67, 2.50, 2.33, 2.17, 3.00, 2.80, 2.60, 2.40, 2.20, 3.00, 2.76, 2.51, 2.28}; // aoz
	vector<double> a  = {0.07, 0.06, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07}; // aoz
	vector<double> y0 = {1.73, 1.73, 1.73, 1.73, 1.27, 1.27, 1.27, 1.27, 1.27, 0.93, 0.93, 0.93, 0.93}; // dE1
	vector<double> b  = {0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17}; // dE1
  // for(auto &x : x0) x -= 0.01;
	// for C16U and C16R //
	double daoz = fabs(x0[0] - x0[1]), ddE1 = fabs(y0[0] - 1.27);
	// C16R
	x0.push_back(x0[0] + daoz); y0.push_back(y0[0]);
	a.push_back(a[0]); b.push_back(b[0]);
	// C16U
	x0.push_back(x0[0]); y0.push_back(y0[0] + ddE1);
	a.push_back(a[0]); b.push_back(b[0]);

	const int nCut = cutls.size();
	for(int I = 0; I < nCut; I++){
		for(int i = 0; i <= n; i++){
			const double theta = 2.*PI*i/n;
			const double x = x0[I] + a[I]*cos(theta);
			const double y = y0[I] + b[I]*sin(theta);
			TCutG *c = cutls[I];
			c->SetPoint(c->GetN(), x, y);
		} // end for over i
	} // end for over I

	for(TCutG *c : cutls) c->Write("", TObject::kOverwrite);

  // Ori cut -- for the full range of the unreacted nucleus //
  const int np = 5;
  double x[np]{}, y[np]{};
  x[0]=2.55913; y[0]=4.16875;
  x[1]=2.55952; y[1]=1.49805;
  x[2]=3.69345; y[2]=1.49805;
  x[3]=3.69345; y[3]=4.17383;
  x[4]=2.55913; y[4]=4.16875;
  TCutG *ori = new TCutG("C16Ori", np, x, y);
  ori->SetTitle("C16Ori");
  ori->SetVarY("pid.dE1()");
  ori->SetVarX("pid.aoz()");
  ori->SetFillStyle(100);
  ori->SetLineColor(kMagenta);
  ori->SetLineWidth(2);
  ori->Write("", TObject::kOverwrite);
} // end main function
